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Abstract 

Excess contributions to the free energy due to interfaces occur for 
many problems encountered in the statistical physics of condensed mat- 
ter when coexistence between different phases is possible (e.g. wetting 
phenomena, nucleation, crystal growth, etc.). This article reviews two 
methods to estimate both interfacial free energies and line tensions by 
Monte Carlo simulations of simple models, (e.g. the Ising model, a sym- 
metrical binary Lennard- Jones fluid exhibiting a miscibility gap, and a 
simple Lennard- Jones fluid). One method is based on thermodynamic in- 
tegration. This method is useful to study flat and inclined interfaces for 
Ising lattices, allowing also the estimation of line tensions of three-phase 
contact lines, when the interfaces meet walls (where "surface fields" may 
act). A generalization to off-lattice systems is described as well. The 
second method is based on the sampling of the order parameter distri- 
bution of the system throughout the two-phase coexistence region of the 
model. Both the interface free energies of fiat interfaces and of (spherical 
or cylindrical) droplets (or bubbles) can be estimated, including also sys- 
tems with walls, where sphere-cap shaped wall-attached droplets occur. 
The curvature-dependence of the interfacial free energy is discussed, and 
estimates for the line tensions are compared to results from the thermo- 
dynamic integration method. Basic limitations of all these methods are 
critically discussed, and an outlook on other approaches is given. 



1 Introduction and Background 

The statistical mechanics of interfaces between coexisting phases dates back to 
van der Waals [T] and Gibbs Some results, like Young's equation '3] for the 
contact angle under which the vapor-liquid interface of a sessile (macroscopic) 
droplet meets the supporting wall in thermal equilibrium, are known since more 
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than 200 years. Understanding the statistical mechanics of interfacial phenom- 
ena is relevant for important problems such as nucleation [il [51 171 151 [TUl ITT] . 
wetting and spreading of fluid films [12lll3l|ll[l5l[T6l[l7l[l8l[l9ll2Ql[ai[22l 
[231 121], metallurgy and crystal growth [5^1 HZ] and diverse technological 
applications that these phenomena have (ranging from oil recovery [55] over the 
efficient deposition of pesticides on plant leaves [5S] and microfluidic devices [5U] 
to the prediction of rain clouds resulting from nucleation of ice crystals in the 
atmosphere [3T], etc.). 

In view of this long history and widespread applications and broad signif- 
icance of this subject it is surprising that still many basic aspects are not yet 
well understood. E.g., regarding interfacial profiles and widths it is still prob- 
lematic to disentangle the "intrinsic" profile and width [TB] from the broadening 
[32l [33l [34] caused by capillary waves [35l |36l [37] . While experimental methods 
exist to measure the interfacial free energy of (macroscopically flat) fluid-fluid 
interfaces and liquid-vapor interfaces, the measurement of interfacial tensions 
involving a solid-phase is difficult [35]: thus such quantities often are extracted 
from contact angle measurements [38j assuming the validity of Young's equa- 
tion [3] (rather than checking it). Since accurate contact angle measurements 
are severely hampered by contact angle hysteresis between advancing and re- 
ceding three-phase contact lines [T^l [HI [131 [23] it is difficult to obtain precise 
experimental information on either the effects of curvature on the interfacial 
free energy of a (nanoscopically small) droplet or on the excess free energy due 
to the three-phase contact line, the line tension [33 HOI EI]- As reviewed by 
Mugele et al. [41] , early attempts to measure the line tensions yielded estimates 
that presumably are several orders of magnitude too large. And although the 
concept of a line tension was already introduced by Gibbs [5] and since then 
has been thoroughly studied (see e.g. Indekeu [32] and Dietrich et al. [33] [33]), 
there are still severe conceptual problems [3S], and there is an ongoing debate 
[46l [47] on the physical significance of the line tension for describing the con- 
tact angle of nanodroplets [48] . Likewise, there has been a longstanding debate 
(see Schrader et al. [33] and Block et al. [SD] for references) on the sign and 
magnitude of the "Tolman length" S [51] introduced to describe the curvature 
dependent surface tension 7(i?), R being the radius of a spherical droplet 

7(i?)=7(oo)/(l + 2^/i?). (1) 

While Tolman [51] originally suggested that (5 is a positive constant and a length 
of molecular dimensions, it soon was recognized that Eq.[l]cannot hold for model 
systems that exhibit a symmetry between the coexisting phases, such as the 
lattice gas model which has particle-hole symmetry, or a strictly symmetrical 
binary (A,B) mixture [52 . Since turning particles into holes (and vice versa) 
turns a droplet surrounded by vapor into a bubble surrounded by liquid, it thus 
changes the sign of R. The requirement that j{R) is invariant against such 
a transformation excludes Eq. [l] for such models, and one finds instead [50] 
another length £ describing a 1/i?^ dependence, 

7(i?)=7(oo)/[l + 2(^/i?)2]. (2) 
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On the basis of a simulation study, ten Wolde and Frenkel [S3] suggested that 
Eq. [2] in fact also holds for asymmetric fluids, such as the ordinary Lennard- 
Jones model for a vapor to liquid transition. However, a suggestion that S = 
for all fluids is at odds with many density functional calculations (see e.g. 
Talanquer and Oxtoby [S3] and Granasy [5S]) that actually 5 is nonzero in 
general but strongly i?-dependent (and 6{R — 00) is slightly negative for a 
droplet surrounded by vapor). Of course, density functional theory is not exact, 
e.g. it ignores the effects due to capillary-wave fluctuations of the interfaces and 
the underlying assumptions of a mean-field character [56] cause the fact that 
the radial density profile of a droplet has a singular behavior (the radius of a 
metastable droplet diverges, although its free energy vanishes, when the state 
of the vapor reaches the spinodal density) , which is spurious for systems with 
short range forces [T] 13 [11] . In view of these problems it is desirable to have 
simulation approaches which do not suffer from these problems. In the present 
paper, we shall review work |50| which suggests that for fluids lacking particular 
symmetries actually a combination of Eqs. [T] and [2] is a reasonable description, 

7(i?)=7(oo)/[l + 25/i? + 2(^/i?)2], (3) 

where the length £ is of molecular dimensions, while S actually is one order of 
magnitude smaller, and negative (for a droplet). This result [50] agrees with 
recent simulation approaches applying rather different methods [57] |5H] and is 
essentially equivalent to a relation d{R) — 6{oo) + const/ R proposed by the 
density functional theory [50] . 

However, also for computer simulation approaches the accurate prediction 
of interfacial free energies has been a challenging problem. One difficulty is that 
one needs to study systems with rather large linear dimensions, to avoid (or at 
least control) finite size effects, that are due to the constraints that boundary 
conditions have on long range interfacial fiuctuations, such as capillary waves 
at interfaces between fluid phases. Even more difficult is the study of solid- 
liquid interfaces [5aiMll6ll|6llMl[M[MllM[6ZllMl[Ml[2Ql[ZDI21 
[75l [76] [77] , since the periodic boundary conditions that are normally applied 
in the directions parallel to the interface may create a misfit resulting in an 
elastic distortion of the crystal structure of the solid, which would lead to severe 
systematic errors. Even when the finite size effects [751|751|771 [751 [751 [50 ] [HT ] 
1551 IMl IMl Uni HZ] are under control, one must consider that interfaces 
are mesoscopic "objects" and their fiuctuations are rather slow, and hence a 
substantial effort of statistical sampling is required. Thus, the present review 
can only be a progress report, and the selection of material is clearly biased by 
the expertise and interests of the authors. 

In the next section, we focus on the estimation of interface free energies from 
thermodynamic integration, using the energy excess of a system with interfaces 
relative to a system without interfaces as observable that is sampled. This 
approach is particularly popular and simple for lattice models such as the Ising 
model [SHI HH HO] , where standard relations for fiuid systems (e.g. based on the 
anisotropy of the pressure tensor [91| ) are not applicable. As an example, we 
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shall discuss the estimation of fint , the interfacial free energy of a domain wall 
between oppositely oriented spins in the Ising model, considering an inclination 
of the interface plane by an angle ■& relative to a (100) lattice plane in the simple 
cubic model [HI]. The method can straightforwardly be extended to obtain the 
excess free energy of the Ising models due to walls [HH IMl (where "surface 
fields" [SSI inni inn HH] may act). In this way, the difference in wall free energies 
of coexisting phases that enter the Young [S| equation, can be directly estimated 
[94l [99] , and this approach can be carried over to off-lattice models of fluids as 
well |1001ll01j . We shall describe such applications in Sec. 3. 

Then wc turn (Sec. 4) to an alternative method, where one samples the 
distribution of the order parameter of a system under conditions where two 
phases can coexist in thermal equilibrium, traversing the coexistence region 
using umbrella sampling techniques |102l 11031 1104] or related methods |105l 
nm \mn mm mm fnOl Iml im Iml HH]. From such studies the Interfaclal 
free energy of flat interfaces of fluid systems can be studied conveniently even 
when it is very small, e.g. in the region near the critical point [1061 11151 11161 
[milllHlinnillinilllllini]. Recently this method has been extended to study 
the free energy of spherical [HI ED] and cyhndrical droplets (and bubbles) [SD])- 

Also an extension to wall-attached droplets has been proposed, which allows 
to test |911 [51] concepts on heterogeneous nucleation jl23) and to provide esti- 
mates for the line tension (Sec. 5). Finally we give in Sec. 6 a brief comparison 
to other methods and an outlook on open problems. 

2 Interfacial free energies of fiat interfaces from 
thermodynamic integration 

2.1 Comments of the preparation of systems containing 
interfaces and the problem of finite size eff'ects 

As is well known, for standard Monte Carlo importance sampling [ 11011113111141 
I124j the bulk free energy _F of a model system is not straightforwardly accessible 
as an "output" of the computation, while the internal energy E is directly 
accessible as the thermal average of the Hamiltonian H, E(T) = (7^)nvTi in 
the NVT ensemble of a fluid (or a lattice system such as the Ising model, for 
instance, where the number N of spins in the lattice and the volume V are 
strictly proportional, rather than being two separate independent variables). 
Using then the standard thermodynamic relation (/3 = l/ksT) 




(4) 
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free energy differences between two states at different temperatures are readily 
obtained from 

l/fcflT 

F{T)^F{To)+ J dp'E{T'), p'^l/keT'. (5) 

l/kBTa 

This simple method has been occasionally used for Ising svstems [1241ll25j noting 
that the free energy is trivially known both for T — > oo, F{T — >■ oo) = E{T) — 
TS{T) -NksT lii2, and for T ->■ 0, F{T ^ 0) = E{Q). In practice, then 
one does not use Tq = in Eq. [s] but rather chooses a nonzero temperature 
To for which the entropy of the system is already negligibly small, so that 
F{Tq) k, E{To) is a sufficiently accurate approximation |124[I125] . As a caveat, 
we note that for off-lattice systems the use of reference states for which the free 
energy is trivially known is clearly more subtle |113| . A general disadvantage 
of Eq. [5] is that for an accurate numerical integration the discretization of the 
temperature interval from Tq to T needs a very large number of intermediate 
states T' , e.g., in the study of surface- induced ordering for the face-centered 
cubic Ising antiferromagnet a number of a few hundred temperatures was used 
|125j , in order to locate the first-order phase transition with a relative accuracy 
of about 3 • 10-5. 

In order to extend this approach to obtain interfacial free energies, one needs 
to choose suitable boundary conditions, in order to prepare a system at phase 
coexistence, with one (or two) interfaces separating the coexisting phases {e.g., 
domains with positive (+) and negative (-) spontaneous magnetization in an 
Ising ferromagnet, see Fig{lJ}-. 

The choice of a negative surface field Hi in the first lattice row (or plane) 
and a positive surface field (of the same absolute strength) Hd = —Hi in the 
last row (or plane), has the disadvantage that the local magnetization m„ in 
the n'th plane or row {n = 1,2, . . . , D) of the lattice will somewhat deviate 
from its bulk value near these free surfaces, over a distance of the order of 
the correlation length ^ in the system [Ml 113 [HS]. The same holds for 
the related "fixed spin" boundary condition (instead of surface fields one has 
spins Si = —1 at all sites in the layer n = adjacent to the first layer on 
the left and spins = -1-1 at all sites in the layer n = D + 1 adjacent to 
the last layer on the right). In fact, for square (or simple cubic) lattices and 
nearest neighbor exchange J the fixed spin boundary condition corresponds 
simply to a specific choice of the surface field {Hi = — J) for the surface field 
boundary condition. No disturbance at the boundaries of the lattice is obtained 
for the antiperiodic boundary condition Fig. [TJd, of course. This choice is only 
possible for systems with a trivial symmetry between the coexisting phases, 
such as the spin-reversal symmetry of the Ising ferromagnet, or the A B 
exchange symmetry of a symmetric binary mixture (A,B). The choices of wall 
fields stabilizing the coexisting phases (Fig. [T^) or the fully periodic boundary 
conditions (Fig. [ij;) are also possible for systems lacking such a symmetry, of 
course. E.g., the choice of Fig. [T^ has been made for the study of interfaces in 
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-periodic boundary condition 
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antiperiodic boundary condition 

-periodic boundary condition 




periodic boundary condition 



Figure 1: Schematic sketch of three possible simulation geometries to study 
interfaces in Ising systems, choosing L x D lattices (in d = 2) or L x L x £) 
lattices (in d = 3), respectively : (a) the "surface field" boundary condition; (b) 
the antiperiodic boundary condition; (c) the fully periodic boundary condition. 
Note that in x (and y, in d = 3 dimensions) periodic boundary conditions are 
applied throughout, to ensure translational invariance in the directions parallel 
to the interface. The interfaces between coexisting domains at positive (-1-) and 
negative (-) magnetization are shown schematically as dash-dotted lines. Note 
that the state (c) is stable if the total magnetization in the system is conserved 
(and chosen zero or close to zero). 
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an off-lattice model for a demixed compressible binary (A,B) alloy where the 
A-rich and B-rich phases have different lattice parameters but the same crystal 
structure [126j . Similarly, the choice of Fig. has been made for the study of 
interfaces in a model for colloid-polymer mixtures |57| . 

In order to extend the application of Eqs. |4] [5] to interfacial free energies, 
one needs to single out the interfacial contribution (which is proportional to 
L or L^, in c? = 2 or d = 3 dimensions) from the bulk free energy (which is 
proportional to LD or Li^D, respectively). This is done by considering suitable 
differences between the energy of a system containing an interface (Fig.[I^,b) or 
two interfaces (Fig. [l];) and a corresponding system without any interfaces but 
the same linear dimensions. E.g., we simulate a system with boundary fields 
Hi = Hd, denoting its energy as £'++(T), while we denote the energy of the 

system of Fig. [TJi as E \-{T). Likewise, we denote the energy of the system 

in Fig. as Eap{T) while a corresponding system with periodic boundary 
conditions throughout (and uniform magnetization, which may be either the 
positive or the negative spontaneous magnetization) has energy Ep{T). Then 
the interfacial free energy -FintCT) is obtained from 

i^int(T)=Fi„t(To)-f j d(i'[E_+{T)-E++{T)], (6) 

l/feflTo 

if the geometry of Fig. [l^ is used, or 

Fint =Fi„t(To)+ j dl3'[EAp{T)^Ep{T)]. (7) 

l/fceTo 

For an Ising ferromagnet we have for Tq — > and an interface perpendicular to 
a lattice direction of the square {d ~ 2) or simple cubic [d — 3) lattice, 

S_+(To) - E++{To) = Eap{To) ~ Ep{To) = ^^int(To) = 2JL''-\ (8) 

due to the "broken bonds" (connecting oppositely oriented spins) across the 
interface. 

Assuming that the two interfaces in Fig. [TJ: are not interacting, a similar 
reasoning (considering the energy difference between the system of Fig. [ij: and 
a corresponding system with also periodic boundary conditions but in a pure 
phase) would simply yield 2i^i„t(T). 

We emphasize that this straightforward approach relies on two basic ingre- 
dients: 

(i) there is a symmetry between the two coexisting phases (denoted by "+" 

and "-" in Fig. [ij, so when we consider differences such as E ^(T) — 

i?+_i_(r), the bulk contribution to the free energy cancels out exactly. 
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(ii) We can reach a reference state (at a temperature Tq) where entropic contri- 
butions are fully negligible, without passing through an intervening first- 
order phase transition when we carry out the integration from l/fcsTo to 

These prerequisites are not fulfilled in many cases of interest, of course, such 
as for off-lattice models. E.g., for a binary symmetric Lennard- Jones fluid (as 
studied by Das et al. [501 [1001 [W] ) Eqs. ||) [f] still would be applicable, but no 
reference state, where the interfacial free energy is known, exists. In fact, at 
low temperatures the system does not stay fluid but undergoes a (first order) 
transition to the crystal phase. Note also, that even for the crystal (due to 
the use of classical statistical mechanics) care is needed in the discussion of 
its entropy at low temperatures, and often one connects via thermodynamic 
integration to an Einstein crystal jl27j . whose statistical mechanics is trivial 
to evaluate. Furthermore, most basic models (e.g., the Lennard- Jones model 
of a simple fluid, or the Asakura-Oosawa model |128j for a colloid polymer 
mixture, etc.) lack a symmetry between coexisting phases. Then no analog of 
antiperiodic boundary conditions exists, while one can still use boundary fields 
that energetically prefer one of the coexisting phases (e.g., this is done in the 
study of interface localization transition using the Asakura-Oosawa model for 
colloid-polymer mixtures |129j ). 

While the simulation geometries shown schematically in Fig. [T] allow the 
estimation of interfacial free energies via thermodynamic integration only in 
exceptional cases, they are widely used to study other interfacial properties, 
e.g. the interfacial profile and its width [13II^,ISH1I11|701|7I1[71[731[71[751 
[761 [IZI [781 [73 [801 [111 [821 [831 (Ml [Ml [Sg IEQ]. However, it is not 

always recognized that these properties do depend strongly on the chosen sim- 
ulation geometry, in particular when the interface is rough (interfaces between 
coexisting fiuid phases always are rough, and the same holds for interfaces in 
the d = 2 Ising model; in contrast, the interface in the d = 3 Ising model is 
rough only for temperatures exceeding the roughening transition temperature 
Tr [131j ). In the geometry of Fig. [l^, the average position of the interface is at 
z = I?/2, but the interface becomes delocalized in the limit D — > cx), similar to 
the case of Figs.[T]D,c, where the z-coordinate(s) of the interface(s) are not fixed. 
This translational degree of freedom of the interface as a whole along the z-axis 
means that there is a logarithmic contribution —ksThiD to Fiat{T) when one 
uses the geometry of Fig. [i]d, and hence there is a logarithmic correction to the 
interfacial tension, for L — > 00, D — )■ 00, 

h,t{T,L,D)^F,^t{T)/{kBTL)^ f,,AT)-\nD/L, d^2, (9) 

fintiT,L,D) = Fi^tiT)/{kBTL^) = /i„t(r) - InD/L^ d - 3 , (10) 

where higher order terms (such as const /i or const/L^, respectively) have been 
ignored. However, for quantities such as the mean square interfacial widths 
w'^ of the interface the size effects are much more dramatic: Fig. [2] gives an 
example j82j . where in the geometry of Fig. [TK the order parameter profile m(z) 
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of an interface between coexisting A-rich and B-rich phases of a model polymer 
mixture was studied |82j . Here m{z) is defined by 



m{z) = [pAiz) - pB{z)]/[pAiz) + Pb{z)] , (11) 

Pa{z), pb{z) being the average densities of A(B) monomers. It turned out that 
m{z) can be described very well by a tanh profile 

m{z) ^mbta.nh[{z - D/2)/w], (12) 

where rrib is close to unity; although Eq. [12] is a standard mean field result 
for "intrinsic" interfacial profiles, where the "intrinsic width" (wq) is a well- 
defined molecular length, Fig. [2] shows that the width w found in the simulation 
has nothing to do with wq: there is a dependence of w on both L and D, 
and on the type of the statistical ensemble used. In the scmi-grand-canonical 
ensemble, the z-coordinate of the interface position may fluctuate around its 
average value zq =^ D/2, while in the canonical ensemble this fluctuation is 
suppressed. Taking (in the semi-grand-canonical ensemble) the limit Z? — >■ cxd at 
fixed L, one hence expects that w'^ cx D^, the interface as a whole may essentially 
fluctuate freely, apart from the immediate vicinity of the confining boundaries 
at z = and z = D. If one takes the opposite limit, L — > oo at fixed I?, a 
phenomenological theory [5T1 rather yields w'^ oc D, while in the canonical 
ensemble w'^ at fixed L reaches finite plateaus when £> — !• oo, cx InL. Both 
laws w'^ (X D, cx In L can be traced back to capillary wave-type fluctuations 
[35J |3S1 [37] . In fact, if one combines capillary wave theory and the "intrinsic 
profile" via a convolution approximation |132| one obtains {D — > oo) 

2 2 , ksT In L fcBrin(g„iax/27r) 
"^"°+4WlT+ 4/.„,(T) ' ^''^ 



where Qmax is a short wavelength cutoff for the capillary wave spectrum. Eq. 13 
shows that from a plot of w'^ vs. InL the interfacial free energy fint{T) can be 
estimated. However, it is not possible to also estimate wg, since in general gmax 
is not known (or perhaps not even precisely defined) . But it has been shown for 



some models [751 [ZSl [EH [Ml [57] that the use of Eq. 13 yields estimates for the 



interfacial tension that are compatible with estimates from other methods. 



2.2 Tilted interfaces and the concept of interfacial stiff- 
ness 

Returning to the Ising model, we now discuss the fact that in lattice systems 
the interfacial free energy does not only depend on temperature, but also on 
the orientation of the interface relative to the lattice directions. Fig. [3] shows 
a generalization of the antiperiodic boundary condition (Fig. ^p) to impose a 
tilted interface. This is done by combining the antiperiodic boundary condition 
(APBC) along the z-direction with a screw-periodic boundary condition (SPBC) 
in x-direction. I.e., the standard periodic boundary condition is not applied at 
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Figure 2: Plot of the mean square interfacial width of an interface between 
coexisting A-rich and B-rich phases of a symmetrical polymer mixture vs. the 
perpendicular linear dimension D for several choices of the parallel linear di- 
mension L (L ~ 32,64,128,256 and 512 lattice spacings, respectively), using the 
geometry of Fig. [TJi. Monte Carlo simulation results for the bond fluctuation 
model of two types of polymer chains with chain lengths Na — Nb — N — 32 are 
shown, at a volume fraction of cj) = 0.5 occupied lattice sites, and a temperature 
T = 0.48rc- Two choices of statistical ensembles are included: the canonical en- 
semble, including identity exchanges in the Monte Carlo move (open symbols), 
marked by (c) and the semi-grand-canonical ensemble (filled symbols), marked 
by (sg). Statistical errors are typically of the order of the size of thy symbols, 
and curves are drawn as guides to the eye only. From Werner et al. [52] . 
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Figure 3: Boundary conditions used for the study of an Lx Lx L Ising system, 
imposing a tilted interface. Antiperiodic boundary conditions (APBC) used in 
the z-direction, an ordinary periodic boundary condition (PEC) is used in the 
y-direction while a screw periodic boundary condition (SPEC) is used in the 
x-direction, involving a shift in the z-direction by Ng lattice spacings. 
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Figure 4: Anisotropic interfacial free energy t{9, T* , L) of the three-dimensional 
simple cubic Ising ferromagnet plotted versus the angle 9, for several values of 
the reduced temperature T* = T/Tc, as indicated; for 6* = 30 two choices 
(L — 64 and L = 32) of the linear dimension of the L x L x L lattice are 
included, to show that at the chosen temperatures finite size effects in fact are 
negligible. Note that r is measured in units of J. From Mon et al. |89j . 



planes z = const, but along the z-axis there occurs a shift by Ng lattice units, 
causing a tilt angle 9 defined via 9 — arctan{Ng / L), while in y-direction still a 
standard periodic boundary condition is used. This combination of boundary 
conditions enforces (at temperatures where the Ising system is in the ordered 
phase) an interface in the system which is inclined (relative to the xy-plane) by 
an angle 9, with the y-axis as the rotation axis. 

Of course, the most general case would involve a SPBC in the y-direction as 
well, involving both a shift Ng_^ for the SPBC in x-direction and a shift Nq^ for 
the SPBC in y-direction. In fact, such a complete study of the angle-dependent 
interfacial free energy Fint{T,9.j;,9y) is mandatory for an understanding of the 
equilibrium shape of a (macroscopic) domain of down-spins surrounded by a 
"sea" of up-spins, in terms of the Wulff |133j construction. However, we are not 
aware of any attempts towards such a complete treatment, and also the work 
of Mon et al. [HS] only deals with the estimation of t{9,T*,L) = Fint{T,9x = 
9,9y = 0)/L2 at low temperatures (T* = T/T^ with [T31| J/ksTc = 0.221655), 
< T* < 0.6 (Fig. |4]). However, for T > Tr (with [HS] ~ 0.54 ± 0.02) the 
angular dependence of t{9, T* , L) is very small [S5^, while for T — > it is quite 
appreciable (Fig.|4|. 
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We use here the notation t{0, T*,L) to alert the reader towards the fact that 
one needs to watch out carefully on finite size effects (although one does not 
detect any in Fig. |4) finite size effects are appreciable near at small 9). In 
fact, the behavior of t{6, T* , L) near is rather subtle, since, for 9 Q 

T((?,r*,L^c^) = T(O,r*)+T'(O,T*)|0|+r"(O,T*)^, r<T^, (14) 

t{9,T*,L^ oo) = t{0,T*) + t"{0,T*){9^/2), T > T*j^ . (15) 

For T < T'^, the interfacial tension has a cusp-shaped singularity at — 0, while 
for T > it is analytic. The derivative t'{0,T*) has the physical meaning of 
a step free energy 

/stop(r) - idT{9, T*,L^ ^)/d9)T'fi=o, (16) 

and the roughening transition is characterized by the vanishing of the step free 
energy (just as the transition in the bulk is characterized by the vanishing of 
the standard interfacial tension between the coexisting phases which become 
indistinguishable at Tc |Slj). The roughening transition is also characterized by 
a diverging correlation length ^r{T) of "height fluctuations" (i.e., "excursions" 
of the local interface height z{x,y) relative to its average z = Zq) |131| . 

^r{T) cx l//step(T) cx exp[c/(Tfl - (17) 

where c is a constant. 

In a finite system the singular behavior for T — >■ Tr is rounded off as soon 
as £,B.{T) becomes of the order of L. We now define for T > Tr the interfacial 
stiffness as p^[T82l[T85] 

k{T*) = [r(0,r)+r"(0,T)]/fcBT. (18) 

In systems with anisotropic but rough interfaces it is k{T) rather than fint/ksT 
that describes the broadening of interfacial profiles by capillary waves {Eq.[l3|}-. 
Fig. [5] compares the temperature dependence of fint{T) and of k{T) for both 
d—2 and d=3 dimensions. It has already been mentioned that (one-dimensional) 
interfaces are rough at all nonzero temperatures, although at T = the inter- 
face (along a lattice direction of the square lattice) clearly is non-rough. Thus 
Tr — can be considered as the roughening transition temperature, and in fact 
ksTniT 0) diverges (Fig.[5|i), while fcsT/int(6' = 0,T ^ = 2), as expected 
from Eq. [8) For interfaces in standard fluid systems, of course, the isotropy of 
space dictates that the interface free energy is completely independent of in- 
terface orientation, and no distinction between the interfacial stiffness (which 
enters as a factor setting the scale for the capillary wave Hamiltonian) and the 
interfacial free energy exists. 
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Figure 5: Schematic temperature variation of the interfacial stiffness ksTn 
and interfacial free energy, for an interface oriented perpendicularly to a lattice 
direction of a square (a) and simple cubic (b) lattice, respectively. While for 
d = 2 the interface is rough for all nonzero temperatures, in = 3 it is rough 
only for temperatures T exceeding the roughening transition temperature Tr. 
For T <Tf{ there exists a nonzero free energy of surface steps (denoted as ksTs 
in this figure) which vanishes at Tr with an essential singularity. While k is 
infinite throughout the non-rough phase k reaches a universal value as T ^ T^. 
Note that k and /int to leading order become identical for T T~, both in 
d = 2 and in d = 3. Note that Tr = in d = 2, and here fcsT/int(T, 6 = 0)/L'^-'^ 
throughout. 
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2.3 Surface excess free energies due to external walls 

The thermodynamic integration approach of Sec. 2.1 can be straightforwardly 
carried over to compute surface excess free energies due to external boundaries 
of a system, such as the waUs of a container holding a fluid, or the free surface of 
an Ising ferromagnet. The Ising case actually is the most instructive and simple 
example, and hence shall be considered first. 

Keeping in mind the interpretation of the Ising model as a lattice gas system 
{Si = —1 and Si = +1 representing the cases of empty and occupied lattice 
sites i, respectively), it is natural to consider thin Ising films with free surfaces 
for which in the first layer (n = 1) a surface field Hi and in the last layer 
{n — D) a surface field Hd act. These surface fields translate in the lattice 
gas interpretation, into short-range potentials attracting the particles of the 
fluid to the walls (which represent the missing layer at n = and n — D + 
respectively) . 

The generic Hamihonian that one uses hence is P l IM l [M l [M l IM l [Wl [55 1 [55] 

H = — J y ' SiSj — Js y ] SiSj — Si 

l^i^j-^ both from the i 

surtaccs" = l °r n = D 

-Hi S,~Hd ' S,^±l. (19) 

i^n—l i^n—D 

Here we consider the nearest neighbor Ising ferromagnet on the simple cubic 
lattice, choosing a Lx Lx D geometry with two free surface layers of size Lx L 
at layers n = 1 and n = D (taking again the lattice spacing as our unit of 
length). Periodic boundary conditions are applied only in x and y directions, 
while there occur no interactions with spins in layers n ~ Q and n = D + \. It 
then is useful to allow in the layers n — 1 and n — D for exchange constants 
Js different from the bulk exchange J. It is well known that the model with 
Js = J (and T > Tr [92]) exhibits at some critical field Hic{T) a critical wetting 
transition |1361I137] . However, choosing Js > J one finds a regime {above a line 
Jsc{T) I where the order of the wetting transition is 1st order rather than 2nd 
order |136l I137j . Since in real systems 1st order wetting transitions are rather 
common [T51 [T71 [THl [T51 [21] and critical wetting is the rare exception [T5] , it is 
of interest to have a model at one's disposal where by variation of a parameter 
{Js/J) one can control the order of the wetting transition. All these wetting 
transitions occur for bulk field _ff = 0, of course. 

In the limit of very thick films {D — >■ oo), where correlations of spins near 
one wall with spins near the other wall can be neglected, the free energy per 
spin /(T, H, Hi,Hd, D) of the model can then be decomposed as [95 ] [96 ] [97 ] [98] 

/(T, H, Hi,Hd,D) = /fc(T, H) + ^fs{T, H, Hi) + ^fs{T, H, Ho). (20) 

Here, fb{T^ H) is the free energy per spin of a bulk Ising system, which depends 
on neither Hi nor Hd, of course. The surface excess free energy of the left 
wall (where Hi acts) is fs{T,H,Hi), while fs{T,H,H]:)) refers to the surface 
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excess free energy of the right wall, where Ho acts. As is well known, wetting 
transitions show up as singularities of the respective surface excess free energies 

We assume now that Hi < and consider the limit _ff 0"'" , so we have (at 
temperatures below the bulk critical temperature T^b) a positive spontaneous 
magnetization rricooxiT) > 0, 

mfc(T, H) = -(9/6(r, H)/dH)T , meocx(T) = mb{T, H ^ 0+). (21) 

We denote the excess free energy of the left wall that belongs to states with 
positive magnetization in the bulk as fs^\T,H,Hi). In the regime of incom- 
plete wetting, fs^\T, H, Hi) then is the excess free energy of a surface where 
the region of positive magnetization in the film extends even close to the left 
wall (although the field Hi would energetically favor a negative magnetization). 
In the wet phase, however, we have a (macroscopically thick) domain of the 
negative magnetization adjacent to the left wall, separated by an interface (as 
studied in Sec. 2.1) from the domain with positive magnetization (which is the 
majority phase in the considered limit H 0+). Consequently, the surface 
excess free energy of a wet surface is 

/r*(T, 0, Hi) - /(-'(T, 0, Hi) + /i„t(T). (22) 

Here /i"^(T, 0, iJi) is the excess free energy of a surface where both the bulk 
magnetization and the surface field Hi are negative, reached for the limit H — )• 
0~ , of course. The wetting transition of this surface occurs when the free energies 
of the wet and the non-wet surfaces are equal, 

/i+)(r,o,i/i) = /i-)(r,o,i/i) + /i„t(T) . (23) 

Note that Eq. [22] explicitly incorporates the fact that the interface separating 
the negative domain (near the wall) from the positive domain is located at 
such a large distance from the wall, that the free energy contributions of the 
wah {fi~\T,0,Hi)} and the interface {/int(r)} can be simply added, without 
considering any interaction forces between the interface and the wall. When 
we consider the Ising-lattice gas analogy, we readily see that Eq. [23] just is 
the standard relation between wall-liquid jwe{T) and wall- vapor ^wv{T) free 
energies at liquid-vapor phase coexistence for a wet wall, 

7^.(7^)= 7-^(7) +7(T), (24) 

where the vapor-liquid interfacial tension is now denoted as 7(T) to make con- 
tact with the notation of Sec. 1. 

For the case of incomplete wetting, the generalization of Eqs. [23] [24] simply 
is the well-known Young |3] equation for the contact angle 0, 

l^,[T)--f^,=-i{T)coBe , (25) 

or 

fi+\T,0,Hi) - /i-)(r,0,Hi) = /int(T)cos0, (26) 
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respectively. The wetting transition is approached from the region of incomplete 
wetting when 9-^0 ^MM\IE\M\niM\MM^^MiM- 

For Ising ferromagnets, spin reversal symmetry implies a symmetry relation 
for the surface excess free energies, 

/i-)(r,0,ffi) = /i+'(T,0,i/i). (27) 

An interesting special case is found for Hi = 0, namely 

/^)(T,0,0)-/i+)(r,0,0) = 0. (28) 

Combining this result with Eg. [26] shows that Hi = corresponds to cos6' = 
0, i.e. 9 = 7r/2. As a result, we conclude that for the estimation of the con- 
tact angle, for which only the difference /i^' (T, 0, Hi) — /i \t, 0, Hi) matters, 
{Eg. 26 28 1- can serve as a convenient initial state for a thermodynamic inte- 
gration. 

In practice one can slightly simplify the problem even further, treating a 
thin film with antisymmetric surface fields, Ho = —Hi, to conclude 

ftHT,0,-Hi) = fi+\T.,0,HD), (29) 

since the surfaces at ti = 1 and n — D can be treated on an equal footing. Now 
one makes use of the relations 



mi = {dfs{T, H, Hi)/dHi)T, mo = -(a/,(T, H, HD)/dHD)T, (30) 

for the desired thermodynamic integrations, since both layer magnetization mi 
and rriD can be straightforwardly sampled. Thus we conclude, using Egs. |27|[30l 

mi 

fi+\T,0,Hi)-fi-\T,0,Hi) = J[mn{T,0,H[)-mi{T,0,H[)]dH[. (31) 



One performs for Eg. [31] a calculation where one studies a thin film with a 
positive magnetization in the bulk and varies the surface fields H'l, H'jj = ~H'i 
from H'l = Q to H'l = Hi in small steps. Fig. |6] shows, as an example, the 
variation of the contact angle 9 ds, a, function of Hi/ J for ksT/J — 4.0 and 
Js/J = 1.0 (a) and Js/J — 1-4 (b). From previous work applying other methods 
|136l 1137] it is known that for Js/J — 1.0 critical wetting at this temperature 
occurs for Hic/J « 0.55 ± 0.015 (while for ksT/ J = 3.0 critical wetting occurs 
for Hic/ J ~ 0.84 ± 0.02) and for Js/J = 1.4 a first order wetting transition 
occurs for Hic/ J ~ 0.3 (due to hysteresis |136) . only a very rough estimation was 
possible). The thermodynamic integration method, using small steps A_ffi/J = 
0.0125 [93 , allows for a much better accuracy, yielded Hi^/ J = 0.318 ± 0.005. 
Near critical wetting, one predicts that [TSl [121 [HI [M] 

A/, = [f^-\T,Q,Hi) - /i+)(T,0,i/i)]//int(+) - 1 (X \HUT) - Hif^" , (32) 
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Figure 6: Contact angle 9 plotted vs. Hi/ J, for the case Js/J ~ 1 (a) and 
Js/ J = 1.4 (b). In case (a), both temperatures ksT/ J = 3.0 (Using a, Lx Lx D 
system with L = D = 60) and fcsT/J = 4.0(L = Z? = 100) are included, while 
case (b) used L = D = 64 and ksT/J = 4.0. Note that at critical wetting 
strong fluctuations and finite size effects render the data very close to Hid J 
unreliable. 
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where v" — 1 according to mean field theory, while for the Ising model v" « 4 
is expected. Since in this region cos^ w 1 — 6*^/2, one concludes further 

B^\HUT)-Hif ,Hi^ HUT) (33) 

as one approaches the 2nd order wetting transition from the incomplete wetting 
regime. From the experience with other methods to study critical wetting with 
Monte Carlo methods ;136.il37j one hence expects that not very close to Hic{T) 
the variation of 6 with Hi is linear, while strong curvature sets in for Hi very 
close to Hic{T). The data of Fig.|6^ are compatible with such an interpretation, 
but do not allow a precise characterization of the singular behavior of 9 near 
Hic{T) due to the limited accuracy of the data used in Fig. |6^ |138j . In the 
case of of first order wetting, we have instead of Eq. [32] 

Af,^HUT)-Hi, e^^\HUT)-Hi\ (34) 

and this relation indeed is borne out nicely by the data. Note that /int(r) in 



Eqs. 26 32 was taken from the work of Hasenbusch and Finn [90 lll39j . see also 



4i. 



3 Estimation of Contact Angles for Off-Lattice 
Models of Fluids 

For off-lattice models of binary (A,B) liquid mixtures, it is also possible to choose 
interactions such that the symmetry against interchange of A and B is strictly 
preserved. This symmetry is equivalent to the spin reversal symmetry of an 
Ising model. The method for obtaining the contact angle of Ising models, that 
was described in the preceding subsection, can then easily be carried over to 
such symmetric models in continuous space. For simplicity, we here summarize 
the approach of a recent "case study" |1001 llOlj that did address a specific 
model, the binary symmetric Lennard- Jones mixture confined in a thin film with 
"antisymmetric" walls, at a temperature far below the critical temperature of 
unmixing for this model. 

In this model, one considers N point particles at positions r^, i = 1, . . . A^, in 
a box of linear dimensions L x L x D, with periodic boundary conditions in x- 
and y-directions, and with impenetrable walls of area LxL each located at z = 
and z = D. The particles interact with pairwise potentials Uaplrij), where a/? 
refers to the type of pair (AA, AB, or BB, respectively), and r^j is the absolute 
value of the distance between the particles, = \ fi —rj\. Starting from the full 
Lennard- Jones (LJ) potential (t>]fj{r) = 4:£ai3[{crai3/rY^ — ((Ta^/r)®], Uai3{rij) is 
chosen identical to previous work addressing the bulk phase diagram |140j . 



(35) 
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while u{rij > Vc) = 0. This choice, Eq. 35 ensures that both the potential 
and the force are continuous at the cutoff-distance r^j = Vc- The potential 
parameters are chosen such that the mixture is fully symmetric, 

CTAA = (^BB = O-AB = CT, ^AA = <^BB = 2eAB = £ , (36) 

and we chose = 2.5a taking units in this section such that a — 1, e = 
1, ks = I- Working at a reduced density p* — pcr^ — Na^/V = 1, the 
critical temperature of unmixing is [140] Tc = 1.4230 ± 0.0005. Considering 
then a temperature T — 1, far below Tc neither the vapor- liquid transition 
nor the liquid-solid transition of the model create problems, the system is a 
dense (almost incompressible) fluid. At this temperature, phase separation into 
coexisting A-rich and B-rich phases is essentially complete (the concentration 
XA = Na/N with N ^ Na + Nb takes a value = 0.97, and x^f^ = 

1 — a;™^°2') = 0.03 by symmetry.) Phase coexistence occurs at chemical potential 
difference Ap = between the particles. 

The "antisymmetric" wall potentials are defined as follows: on the A-particles 
acts a purely repulsive wall at both walls, and in addition an attractive poten- 
tial (of strength e^) only at the left wall (z = 0). On the B-particles, the same 
repulsive potentials acts on both walls, and in addition there is an attractive 
potential (of strength Ca) only at the right wall (z=D). Thus 

where an offset S — (t/2 is used, and = e/15. 

One then simulates the thin film in the semi-grand-canonical ensemble under 
phase coexistence conditions (A/x = 0) and varies e^. The difference in wall free 
energies JwA^JwB in the regime of incomplete wetting (considering, with no loss 
of generality, the A-rich phase in the bulk) is then found by a thermodynamic 
integration method, inspired by the treatment of Sec. 2.3. We write the free 
energy of the thin film in terms of the partition function Z as follows (/3 = 
l/keT) 

F = -kBT\nZ= -ksTln J dX cxp | - (SHbiX) - miiX) - l3eaHl{X) 

(39) 

where X stands for the coordinates of all the particles, Hb{X) describes all the 
particle interactions, Hl^{X) the energy of all the particles due to the repulsive 



part of Eqs. 37 38 and ea'H'^{X) the energy of all the particles due to the 



attractive part of Eqs. 37 |38| Prefactors in Z that are unimportant in the 



present context have been not spelled out. Note that (AT) can be written as 



'H-[X)^L^'^-^)[j PA{z,X){ — fdz + J pB{z,X)i ^^"^_J 'dz] (40) 
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Here pa{z,X) is the density of A-particles in configuration X in the (infinitesi- 
mally thin) shoe L^dz from z to z + dz, and X) is defined similarly. Note 

that due to the repulsive wall potentials, pA^z^X) and pb{z,X) are essentially 
zero outside the interval < z < D, over which the integrations in Eq. [39] are 
performed. 

Since for D — >■ oo, the free energy F can be decomposed as 



{df:''/dea)T,A>^=0 



27rp* 



D 



dz{ — )'{pA{z)). 



(41) 



and 



2-Kp* 



dz{ 



D 



\pb{z)). 



(42) 



Here the notation (. . emphasizes that the averages are sampled at a nonzero 
value of ta in Eqs. [37 38 Eqs. 41 42 are the generalizations of Eq. 30 for the 
Ising (lattice gas) model - if we chose in the Ising model a "surface field" Hif(n) 
which does not only act on the spins in the first layer n = 1, but on all layers, 

by "i„/(n), of course. 



we would have to replace mi in Eq. 
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Thermodynamic integration of Eqs. 

flcft/, ^ _ Aciti, _ n\ , ^""^ 



41] [42| with respect to then yields 



/r(o = /r(ea = o) + 



de'a / dz{ 



(43) 



/r'^^'(ea)=/r«"(ea=0) + 



right/ _ 



27rp* 



D 



de'a dz{- 



p — z 



-f{PB{z)). 



(44) 



Remember that so far we were considering solely wall free energies of the 
A-rich phase. But due to the symmetry of our model the wall free energies of 
A-rich and B-rich phases are related; similar to Eq. [27]we have 



B — rich 

Thus the desired difference can be written 



(45) 



A-rich 



/r(£a)| 



B-rich 



fs (£a)|A-iich ~ fs^ i£a)\A-rich — 
D 

= (^)p< / 


(^)^-(p^(z)).,( / 

2 + " U + — z 



(46) 
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38 plotted vs. 



Figure 7: Contact angle of a symmetrical binary (A,B) Lennard- Jones mixture 
1.0, exposed to wall potentials Eqs. 37 
of the attractive wall potential, 
is predicted to occur at « 0.24. 



{Eqs. [35) [36]1-, at T 

the strength 



ea, analogous to Eq. 
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A fist-order wetting transition 
The data are compatible with a relation 
(broken curve) 



Both statistical averages {pA{z))^t^, {pb{z))^i^ are sampled in the same (A-rich) 
phase only. Of course, — IwB — for = 0, implying once more a contact 
angle of 90° for this "neutral wall" situation. Subtle conceptual problems [IS] 
concerning the proper identification of individual wall free energies ^wA,1wB 
due to the freedom to define exactly where in our off-lattice model the walls are 
located {z — Q and z = D or z = —5 and z = D ^ 6) will however leave this 
difference unaffected. 

If the interfacial tension ^ab between the coexisting A-rich and B-rich phases 
is known independently (see Sec. 4), Young's equation cos 9 = ['^wA — 1wb]/i{T) 
can be used to predict the dependence of the contact angle on (Fig. 7). 
In this model it is also possible to extract estimates of 9 directly from the 
observation of slab configurations |1001I101] equilibrating thin films (for the same 
choice of bulk and wall potentials as described above) in the canonical ensemble, 
using xa = xb = 0.5, and starting with a slab configuration with perpendicular 
walls at xi = L/A and X2 — 3L/4 as an initial condition (the B-rich phase 
is then located in between xi and X2)- Although the necessary equilibration 
of the interfaces is a time-consuming procedure, accurate estimates for 9 at 
various values of Ca < ea,c can be obtained [1001 llOlj . however choices where 
9 is small cannot be studied (due to the need to make both L and D rather 
large). Somewhat surprisingly, it was found that the "macroscopic" contact 
angle (as obtained from the method as described above, which avoids the study 
of phase coexistence in the thin film geometry, focusing on estimating suitable 
observables exclusively in the single-phase region of the system) agrees with 
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the contact angles found from direct observations of phase coexistence on the 
nanoscale very well [1001 llOlj . 

Unfortunately, an extension of the above method to systems lacking this 
particular symmetry does not seem to be straightforward: surface potentials 
typically in a thin film will produce a shift of the parameters where phase 
coexistence occur (by a mechanism analogous to "capillary condensation" ) . 
E. g. for an asymmetric binary mixture the bulk phase coexistence occurs at 
a nontrivial value AficoexiT) and in the thin film (even with a purely repulsive 
potential at both walls) will be shifted to a value A/icoox(T', D). Thus, this case 
also corresponds to a nontrivial (a priori unknown) value of the contact angle 
at both walls, implying that the interface separating both coexisting phases 
and running from the wall at z = to the wall at z — D is curved. Such 
curved interfaces occur also in symmetric mixtures with symmetric rather than 
"antisymmetric" wall potentials (see e. g. |141j for examples), or when one 
considers "liquid bridges" at vapor-liquid phase coexistence in capillaries (e.g. 
|142j ). and are hard to analyze. 

An elegant way to avoid this problem, allowing the direct observation of 
contact angles of flat interfaces between coexisting phases that lack particular 
symmetries, has been proposed by Dimitrov et al. |143j . One chooses a thin 
film geometry where phase coexistence between liquid and vapor at coexistence 
pressure of the bulk is enforced by coupling the system to an external liquid 
reservoir held at this pressure. The wall at z = has two different parts: for 
< X < L/2, it has a wall potential that is strongly attractive to the fluid 
particles (leading to complete wetting of the liquid), while ior L/2 < x < L 
it is repulsive (leading to drying), and a periodic boundary condition is used 
only along the y-direction. The liquid reservoir is at a; < 0, while at a; < L 
there is another (repulsive) wall. The other wall at z — D is homogeneous for 
< a; < L, and there the wall potential is chosen such that incomplete wetting 
conditions occur. In this way, one can stabilize a flat interface, pinned at z = 
along the line x = L/2, and "hitting" the wall at z = £> under a nontrivial 
contact angle (which is the observable of interest) that depends on the potential 
that acts on this wall. Only near z — does one expect some curvature of 
the interface, essentially one creates a flat interface inclined under the contact 
angle, which can be varied over an extended range |143j . 

4 The order parameter distribution and what 
we can learn from it about phase coexistence 

In this section we consider systems that exhibit a critical temperature Tc such 
that for temperatures T < Tc the system can coexist in two phases distinguished 
by different values of a scalar order parameter: e.g., in an anisotropic magnet 
(as described by the Ising model) the order parameter is the magnetization 
per spin to, and domains with magnetization ±TOcoox can coexist (tocoox being 
the absolute value of the spontaneous magnetization); in an (incompressible) 
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binary (A,B) mixture the order parameter is the relative concentration of one 
species, say A, defined as xa = Na/{Na + Nb) where Na, Nb are the particle 
numbers of species A and B. We then assume that for T <Tc a. miscibility gap 
occurs, described by the coexistence of two phases with different concentrations, 

^^'^ ■^A(2)- consider here explicitly again only the case of a strictly 
symmetric Lennard- Jones mixture (details about this model were already given 
in the previous section) for which x'^^°^-^ = l — x'^^°^y Fig. 8 shows, as an example, 
the phase diagram of the model that we study (it has a critical temperature at 
Tc = 1.4230 ± 0.0005, while we study phase coexistence at T = 1.0 where 

~ 0.030, so one has coexistence between almost pure A and pure B). 
Finally, we deal with the simple Lennard- Jones fluid [321 [SO] where the order 
parameter is the particle density p — N/V (V being the volume of the simulation 
box, which we take as a cube of linear dimension L with periodic boundary 
conditions), and for T < vapor (at density pi,)and liquid (at density pi) 
can coexist. The model that is studied is defined by the Lennard- Jones (LJ) 
potential u{r) depending on the distance r between the particles, 

u{r) = 4e{(CT/r)i2 - {a/rf} + C, r< r^, (47) 

u{r) = 0, r> rc, (48) 

choosing units such that the depth of the potential well (e) and its range (a) 
both are unity, e = l,cr = 1, the cutoff rc — 2.2^/^(t, and the constant C is 
chosen as C = 127/16384 so u{r = = 0. For this model = 0.999 [133] 
(also choosing Boltzmann's constant unity). 

While for the Ising model spin reversal symmetry ensures that phase coex- 
istence occurs at magnetic field H = 0, and in the symmetric binary mixture 
the symmetry against interchange of particle labels A and B ensures that phase 
coexistence occurs at chemical potential difference A/j, = /i^ — /is = 0, no such 
symmetry exists for the one-component Lennard- Jones fiuid, and phase coex- 
istence occurs at a nontrivial value /icocx(T') of the chemical potential p, of the 
particles. Carrying out simulations in the grand canonical {pVT) ensemble and 
recording the probability distribution of the density Pl(p, T) one finds Pcocx(T) 
from the equal weight rule {in the region of p where Pl (p, T) exhibits two peaks 
one near Pv{T) and the other near p^(T)}. [1451 1146j . 

Having obtained the distribution functions PLHri'm) for the Ising model 
(at H — Q),Plna^it{xa) for the symmetric binary mixture (at A/z = 0), and 
Pl^it{p) for the LJ fluid (at p = PcocxiT)), we associate with these distributions 
effective free energy functions, defined as follows (the dimensionality d is 3 
throughout in the following; and for simplicity, we use the symbol / throughout 
for these densities of thermodynamic potential, irrespective of the system) 

/L(m,T) = ~ikBT/L'')ln[PLHTim)/PLHT{m,oe.)l (49) 

fUxA, T) = -{kBT/L'') HPLNAM^A)/PLNA^.T{xT'')] , (50) 

and 

/i(p,T) = -{kBT/L") HPlMp)/PlMPv)] + W • (51) 
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Figure 8: Phase diagram of the symmetric binary (A,B) Lennard- Jones mixture 
with p* = 1 in the plane of variables T and relative concentration xa — Na/N 
of the A particles, for fixed total particle number N — 6400. The cross shows 
the critical point, taken from |140| . The horizontal broken line indicates that 
phase coexistence is studied for T = 1.0. From block et al.fSU]. 



Typical examples for these effective free energy functions are shown in Fig. [9j 
One can see a characteristic double-well shape, with a "hump" in between that 
is more or less flat in the center, and depends clearly on the linear dimension of 
the system. 

The double well-shape of these effective free energies may be considered as 
somewhat reminiscent of the double well shaped mean field free energy density 
/*=f^(m, T) obtained when one treats the Ising model in molecular field approx- 
imation {and related mean field free energy densities f^'^{xA,T),f'^^^{p,T) 
for the other systems}. We emphasize, however, that such an analogy would 
be completely misleading: while these mean field free energy densities do not 
depend on L and describe inside the two-phase coexistence region, homogeneous 
states that are interpreted as being metastable (up to the inflection point, the 
so-called spinodal) and unstable (inside of the spinodal) , we deal here with the 
statistical mechanics of finite systems in inhomogeneous states in full thermal 
equilibrium. The whole structure that is seen in these curves inside of the 
two-phase coexistence region all comes from interfacial effects [5]. 

This statement is easily verified when one examines snapshot pictures of 
the typical configurations contributing to the average in various regions of this 
effective free energy density. As an example. Fig. [lO] presents configuration 
snapshots |100) of the symmetric binary {A,B) mixture for xa = 0.15, 0.27, 
and 0.50, respectively. In the snapshot pictures the positions of the A-particles 
are indicated by dots, the positions of the B-particles are not shown at all, for 
the sake of clarity. One clearly recognizes that both droplets of approximately 
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Figure 9: (a) Effective free energy per spin /L(m, T) of Lx Lx L Ising lattices 
plotted vs. m at hgT/J = 3.0 for three choices of L, as indicated, (b) Effective 
free energy fLixA, T) of a symmetric binary Lennard- Jones mixture in a i x L x 
L cube with periodic boundary conditions for several choices of L and T — 1.0. 
From Block at al. [50] . (c) Effective free energy fL{p,T) of a Lennard- Jones 
fluid {Eqs. [48| [ at T = O.TSTc, for six trices of L as indicated. 





(a) 




Figure 10: Configuration snapshots of a symmetric binary Lennard- Jones mix- 
ture at T = 1.0 for L x L x L systems with periodic boundary conditions with 
L = 24 and three different choices of the concentration: xa = 0.15 (a), 0.27 
(b) and 0.5(c). Only A particles are shown by dark dots, B particles are not 
shown for clarity. In case (a), one can recognize an (almost) spherical droplet 
surrounded by supersaturated "vapor" of A particles; case (b) shows an almost 
cylindrical droplet, connected into itself via the periodic boundary condition; 
case (c) shows a slab configuration (the A-rich phase is separated from the phase 
poor in A by almost planar interfaces). 
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spherical and cylindrical shape occur, as well as slab-like configurations bounded 
by two (approximately planar) interfaces. The latter two configurations are 
stabilized by the periodic boundary conditions. Note that for a.n L x L x L box 
the cylinder axis can be oriented along any of the x,y,z coordinate axes, and 
the same symmetry holds for the orientation vector that is perpendicular to the 
interfaces in the slab state. Of course, these spherical, cylindrical and slab-like 
domains exhibit these regular shapes on average only: the sampling that is done 
does not at all constrain the fluctuations of these domains, at every value of xa 
all configurations that are compatible with the chosen concentration contribute 
according to their statistical weight. 

By choosing rectangular boxes one can constrain the orientations of the 
cylinders or slabs, respectively: choosing a volume L x L x D with D < L 
the cylinder axis will be predominantly along the z-direction (since a cylinder 
of radius R and height D then costs a surface free energy proportional to a 
surface area 2RttD instead of 2RttL). Likewise, choosing a volume L x L x D 
with D > L the domain walls in the slab state will predominantly be oriented 
perpendicular to the z-axis, since the surface free energy cost (proportional to 
L^) then is less than for the other two orientations (for which it is proportional 
to LD > L^). 

We now comment on the flatness of the hump in Fig. [9}d in the region when 
the slab configuration prevails. If we assume that the two interfaces that sep- 
arate the A-rich phase from the B-rich background do not interact with each 
other at all, the free energy cost to form them is simply 27^5 (L)L^, where 
^ab{L) is the interfacial free energy density (since the periodic boundary con- 
ditions constrain interfacial fluctuations, e. g. the spectrum of capillary waves is 
discretized, we allow some weak dependence of ")ab{L) on L, requiring however 
that 7_ab(oo) exists). 

The statement that the droplet shapes are (on average) either spherical or 
cylindrical have tacitly implied that the interfacial free energy does not depend 
on the direction of the vector perpendicular to the interface. This is the case for 
off-lattice fluid-fluid interfaces (such as the interfaces of the model binary fluid 
mixture studied in Fig. 10 ) or liquid- vapor interfaces. However, for the Ising 
model, due to the underlying lattice structure, such an isotropy of interfacial 
properties holds only very close to the critical point. In fact, as discussed al- 
ready in Sec. 2.2, for temperatures below the roughening transition temperature 
Tfj of the (100) interface the droplet surface contains strictly flat regions. While 
at T = the shape that minimizes the surface free energy at given interface 
orientation in the Ising model simply is a cube, at nonzero T < the determi- 
nation of the minimum free energy droplet shape is a nontrivial problem [133j , 
since the interface contains both planar parts ("facets") and curved regions, 
and since facets appear only for T < Tr, Tr can also be interpreted as faceting 
transition temperature. Similar complications due to the anisotropy of the in- 
terfacial free energy are also important for crystal-fluid interfaces. However, we 
shall not further consider the problem of anisotropic interfacial free energies in 
this paper. 

Another important problem concerns the transitions between the various 
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Figure 11: Extrapolation of the interface tension jAsiL) of flat interfaces 
between A-rich and B-rich phases of Lennard- Jones mixtures, coexisting in L x 
L X L systems at ksT /e = 1.0 in a slab geometry (as discussed in Figs. |9]jb), 
^c)). The shown straight line fit yields 7as(oo) = 0.722 ± 0.002. From Block 
et al. i50j. 

regimes in Figs.[9|and |10[ The transitions between slabs and cylinders, cylinders 
and spheres, and states with spherical droplets and without it are not sharp (for 
finite L) but rounded, when xa is varied. This observation is consistent with the 
fact, of course, that "sharp" phase transitions (i.e., associated with a singular 
variation of the free energy of the system) can occur in the thermodynamic 
limit only, L — ^ cx). In this limit, the transition (at x\) from the homogeneous 
one-phase state to a state containing a droplet moves towards the coexistence 
curve, x\ ~^ x'a°i^ while the droplet to cylinder transition occurs at a volume 
fraction ^ of the A-rich phase given by 

^'^"P-'^y' = (47r)/81 w 0.155, ^ {x a - x^Xf) / {I - ^x^XT) (52) 
and the cylinder to slab transition occurs at 

^cyl -slab ^ ^cocx ^ ^cyl-slab^^ _ 2xXf), ^"^^-''^^'^ = I/tT « 0.318 . (53) 

Due to the symmetry of the system with respect to the interchange of A and 
B, analogous transitions occur at 1 — \li^'°P"^y' and 1 — vj/'^yi-'^iab^ Yoy finite L, all 
these transitions are rounded, and hence when one analyzes the interfacial free 
energies of spherical or cylindrical interfaces (or even flat interfaces encountered 
in the slab states), it is important that one stays off from the regions where these 
rounded transition occur. Of course, in the strict sense this means that always 
only the approach towards the thermodynamic limit is a well-defined problem, 
both for curved and for flat interfaces. In fact, for a spherical droplet there 
is always a nonzero (albeit for large L in practice negligibly small) probability 
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that the droplet makes a transition to cyhndrical shape or that it evaporates 
completely. This caveat is not just a matter of principle, but actually prevents 
the study of extremely small droplets, that contain just only a few particles 
only, also in practice - such small droplets are not well distinguished from other 
fluctuations that occur in the homogeneous phase. 

Fig. [TT] shows, as an example, the estimation of the interfacial tension jab = 
jAsiL — co) for the binary LJ mixture, using an extrapolation ofjABiL) versus 
1/L. It is seen, that in practice the dependence of jab{L) on L can be rather 
weak, and a very good accuracy can in fact be reached. However, it should 
be noted that the precise nature of finite-size corrections to jab is not fully 
understood. When this method of estimating interfacial tensions from /^(m, T) 
was introduced almost 30 years ago for the Ising model in d = 2 and d ~ 3 



dimensions |115j , an alternative formula including a logarithmic correction was 
also suggested 

/rx ^ \ b\nL 
7(L)=7(oo) + - + — ^, (54) 

with two undetermined constants a,b. However, lacking theoretical guidance on 
the values of these constants, a very large range of L would be required to allow 
for an unambiguous estimation of the three constants 7(00), a and b from an 
empirical fit to simulation data on j{L). 

We note that despite these difficulties, there is ample evidence both in d = 2 
(where one can compare to Onsager's exact solution for 7(00) |143) ) and in d = 3 
that the method illustrated in Fig.[Tl]is practically useful [Ml ESI [Z3 [Ml IMl 
[TTHl [TT71 [TT8l [TT9l [T20I [ml [T22] . In fact, for accurate estimations of 7„^ of 
the vapor-liquid interface tension this method has become the method of choice 

[igiimiTTTifngi iTmfT^ fnniT^. 

Thus, while it early has been recognized that /L(m, T) {as well as fL{xA, T) 
or fL{p,T), respectively} contains valuable information on the interface free 
energy of planar interface in the respective models, only recently it has been 
exploited that one can obtain also the interface free energy for the case of 
curved interfaces. Fig. [T2| explains the basic idea for the case of the fluid binary 
mixture. One first introduces the variable conjugate to the order parameter 
(m, XA, or p, respectively) as a derivative of the free energy with respect to the 
order parameter: 

HLim,T)^[dfLim,T)/dm]T, (55) 
^A^iL{xA,T) = [dfLixA,T)/dxA]T, (56) 
and ^ 

f,L{p,T) = [dfLip,T)/dp]T. (57) 

In the thermodynamic limit, H]^(m,T) converges to an L-independent func- 
tion H{m, T) for m < — TOcocx; this is just the description of the magnetic equa- 
tion of state of the Ising magnet (or the equation of state of the binary mixture 
or fluid) in the bulk one-phase region, of course. Analogous statements refer to 
Apl{xa,T) and pl{p,T). Fig. 13 shows specific examples for these functions, 
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Figure 12: Schematic explanation of the generaUzed lever rule and its use to 
obtain droplet free energies for a symmetrical fltiid binary (A,B) mixttirc. At 
the same chemical potential difference (horizontal brokcm line) three states 
can be identified: a pure B-rich phase at concentration a;^''' > x'-''^°^, a pure A- 
rich phase at concentration > x'^J^f' = 1 — x™'''', and a state containing 
an A-droplet of radius i?, surrounded by B-rich background (which must also 
have the concentration x^^), at average concentration xa- Defining R by the 
condition that the concentration in the droplet equals x^, observation of Ax 
yields information on i?, while observation of A/ yields information on the 
surface free energy of the droplet. 
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Figure 13: (a) Ising model isotherms at ksT/J = 3.0 plotting HL{m,T)/kBT 
versus magnetization m for several choices of L as indicated. From Winter et al. 
[SI], (b) Isotherms of the binary symmetric Lennard- Jones mixture at T = 1, 
yO = 1, plotting A^l{xa,T) /UbT vs. concentration xa- Data refer to a single 
run at each size L, to illustrate the typical noise level (200 Monte Carlo steps per 
particle have been used for each window of the successive umbrella sampling). 
For the final analysis, five such runs wife averaged together. From Block et al. 
[50] . (c) Isotherms of the simple Lennard Jones fluid {Eqs. 48|- at T = O.TSTc, 
plotting pll{p,T) /ksT versus the density p, for six linear dimensions L, as 
indicated. Note that lengths are given in units of cr(= 1). From Block [149]. 



defined in Eqs. [SSpT) as obtained from simulations [i^ [5U1 IMl I149j . The data 



show that the maxima and minima of these curves (which correspond to the 
droplet (bubble) evaporation/condensation transition [1501 llSlj I152j . as noted 
above) indeed become sharper as L increases, and move towards the coexistence 
curve. Arguments have been presented [15011151] that in d = 3 dimensions both 
the distance from coexistence and the height of these extrema scale as L^^/^ 
and some numerical evidence for this relation can be found in Schrader et al. 
|148j . For d = 2, an analogous relation (scaling with L~^/^) has been estab- 
lished rigorously |152j . Thus, when the limit L — >■ c» has been taken, no trace 
of this transition (nor of the transitions from spherical to cylindrical domains 
(Eq. |52| or from cylindrical domains to slabs (Eq. |53| , respectively) is left in 
the isotherms, while for finite L all these transitions are subject to finite size 



rounding, as is obvious from Fig. 13 However, for large L the relative extent of 
the rounding is small (see Ref. [151j for a phenomenological discussion), and 
hence for large L one can identify parts of the rounding shown in Figs. [9] and 
[13] where the droplet-vapor coexistence is not affected by these transitions, and 



hence the analysis outlined in Fig. 12 becomes applicable. Similarly, for volume 
fractions of the minority phase in the region \I/'i''°P-'=yi < ^E* < vp^yi-siab^ 
can use an analysis analogous to Fig. [12] to extract the surface free energy of 
cylindrical interfaces [SD] as long as one does not use volume fractions ^E" close 
to the boundaries of this interval. 

We here indicate only for the case of the binary mixture, how the surface 
free energy of droplets is extracted, inspired by the construction indicated in 



Fig. 12 It is implied that the volumes of the two phases Vi (droplet) and Vn 



(surroundings) and their particle numbers Na.i, Na,ii,Nbi, Nbii are additive 

V ^Vi + Vii,Na^ Nai + Naii , Nb = Nbi + Nbii • (58) 

Thus the interface in Fig. [12] is defined infinitesimally thin (so there is no excess 
volume associated with the interface) and its location is placed such that there 
occurs no excess of the particle number of either component. Note that we 
have assumed a strictly symmetric and practically incompressible mixture, so 
that a description in terms of a single order parameter density xa suffices. 
Unfortunately, the generalization of our procedures to compressible mixtures 
where two order parameter densities p and xa need to be considered is not 
straightforward. 

Fig. [12] then implies, when we assume that the average shape of the droplet 
is a sphere of radius R, that 

Ax = x'X"^ -XA = (47ri?V3L3)(l - 2xT'") ; (59) 

remember that a;^^ is the concentration of the surroundings (supersaturated in 
A particles) of the droplet. Since the free energy density of this phase and the 
bulk term of the free energy density of the droplet are equal. A/ is entirely due 
to the surface free energy of the droplet, and hence, 

Af = AnR'-fAB{R)/L\ (60) 
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Figure 14: (a) Plot of j(oo)/'-f{R) — 1 versus for the Ising model at 

ksT/ J = 4.0. Here 7ab(oo) was taken from the work of Hasenbusch and Pinn 
[3D]. Different symbols show different choices of the linear dimension L, as 
indicated, (b) Plot of ^AB{'x)hAB{R) ~ 1 versus for the symmetric 

binary Lennard- Jones mixture at T = 1, p = 1. Here, ^ab was taken from 
Fig. [TT] Different choices of L are used for different parts of the curve. From 
Block et al. [50]. 
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where we have denoted the interfacial tension of a droplet, which has the 
(equimolar) radius R, as jAsiR), for the case of the symmetric binary mix- 
ture. Of course, Eqs. [59] [60] can be straightforwardly transcribed to the Ising 
model (if the approximation that the interface tension is independent of inter- 
face orientation is accurate). In the case of the vapor-liquid transition, which 
lacks any symmetry between the coexisting phases, of course, droplets of radius 
R and bubbles of radius R may give rise to different surface tensions Jvl{R) 
in both cases. As is well-known, this difference is related to the existence of a 
nonzero Tolman length j51) . as discussed already in the introduction. 

Using this method outlined in Fig. [12] the curvature-dependent interfacial 
tension has been obtained both for the Ising model [IHl IM] , the binary symmet- 
ric Lennard- Jones mixture [SU] and the simple Lennard- Jones fluid [SHI 1148] . 



Fig. 14 presents plots of 7(oo)/7(i?) — 1 (Ising model) and 7AB(oo)/7AB(i?) — 1 
(binary mixture) versus 1 /R^ ; the linear variation demonstrates that these data 
are compatible with Eq. [2] Of course, due to the need to use, for each choice 
of L, only those data in Fig. [13] which are not affected by the droplet evapo- 
ration/condensation transition nor by the sphere to cylinder transition of the 
droplet, each choice of L can yield only a small part of the desired curves 7(i?) 
and 'yAsiR)- Ideally, these small parts should superimpose such that a unique 
master curve results. Actually, there clearly occurs some scatter rather than 
such an ideal perfect superposition, but nevertheless the desired master curve 
is reasonably well defined, yielding 

7(oo)/7(i?) - 1 « 10/i?^ Ising model, (61) 

7AB(oo)/7yis(^) — 1 ~ 2.2/i?'^, binary mixture. (62) 



While for the systems shown in Fig. 14 the existence of a Tolman length [5T] 



is precluded simply due to the symmetry between the coexisting phases, for the 
simple Lennard- Jones fluid there is no such symmetry between vapor and liquid, 
and a Tolman length S should exist. Since the radius of curvature is negative for 
a bubble, the sign of the correction in Eq. [l]for droplets and bubbles must differ. 



Fig. 15 shows plots of "fyi{R)/^yi{oo) — 1 vs. 1/R for two temperatures. Indeed 
one sees that this curvature correction for droplets differs from its counterpart 
for bubbles, ruling out that Eq. [2] holds. On the other hand, in the available 
region the data are not compatible with a simple linear relation in 1 /R {Eq. [ij- 
either, while a flt to Eq. [3] which involves two characteristic lengths S and i is 
possible. Taking these data together with results for cylinders, for which one 
expects instead of Eq. [3] a relation 

7(i?) = 7(f»)/l[l + S/R + 2{£,/R)% (63) 

where the expectation is that the same length 6 in the leading 1/R correction 
must occur, while the lengths for spheres {£) and cylinders {£c) in the quadratic 
corrections may differ. Estimates for S « —0.1a were extracted |50j at both 
temperatures that were studied, while the lengths i, £c were both found to be 
of order cr, as for the binary symmetric LJ mixture. The fact that 6 is nega- 
tive, and its absolute value is only of the order O.lcr, is compatible both with 
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Figure 15: Plots of jye{oo)/jyg{R)~l vs 1/R for spherical droplets and bubbles 
for the LJ fluid at (a) T = 0.78Tc, and (b) T = 0.68Tc. Fits to functional forms 
y — ax + bx^, with adjustable constants a,b as quoted in the figure, are included. 
From Block et al. EOl. 
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Figure 16: Schematic sketch of an isotherm H^ jj versus m in the Ising model 
in the L x L x D geometry where in the first layer (n = 1) a positive surface 
field Hi acts, while in the last layer {n = D) a negative surface field H]j = —Hi 
acts (upper part). The strength of these surface fields is assumed to be in the 
regime of incomplete wetting {Hi < Hi^iT)}, and the temperature is below 
criticality (T < Tc). The lower part shows the associated effective free energy 
fL,D{'m, Hi,T) of the system. For both HL^oi'm, Hi,T) and fL,D{'m, Hi,T) 
one can distinguish different regimes: for m < —mt{m > rrit) pure phases 
of predominantly negative (positive) magnetization occur; for —rrit < m < 
—rrit' (m'f < m < mt) wall attached sphere-cap shaped droplets of positive 
(negative) magnetization occur, while for —m[ < m < — (m" < to < mj) 
the shape of the wall-attached droplets has changed from sphere cap to cylinder 
cap (not shown). In the region — m" < m < to", the walls are planar, and thus 
H]^{m, Hi,T) = in this regime. The phase coexistence in the bulk occurs for 
TO — ±TOcoox- The dotted horizontal line in the upper part indicates that the 
(pure) state with magnetization to' at the ascending part can coexist not only 
with a pure state at the other ascending part (to') but also with a mixed-phase 
state (m) at the descending part of the curve iJi.£i(TO, i/i, T). Note that the 
transition at the values ±TOt,±TO(, and ±m" (highlighted by broken vertical 
lines) are not sharp but exhibit finite-size rounding. 
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Figure 17: (a) Plot of free energy hump /^(m w 0,i?i = 0,T)L/{2kBT) for 
the Ising model at keT/J = 3.0 (a) and of /^(xa « 1/2, = 0,T)L/2kBT for 
the binary Lennard- Jones mixture at 

2^8, 1.0 (b) versus 2/D. Linear dimension 
L — 40 lattice spacings in case (a) and L = 30 Lennard Jones units in case (b). 
Part (a) is taken from [99] , part (b) from [50] . Note that the ordinate intercept 
in case (a) was fixed at the literature value |139[ 7 = 0.434, and in case (b) at 
the estimate ^ab = 0.72 (Fig. [Tl]). The slope of the straight lines yields line 
tension estimates t{0 = tt/2) « -0.26 ± 0.01 (a) and -0.52 ± 0.01 (b). 



density functional calculations [SO] and with estimates extracted from different 
approaches [571 [SS]. However, Anisimov |153j has derived that near the critical 
point one expects S to diverge as 

6 = -S{1 - T/T.f-'' , (64) 

where /3 ~ 0.325 and v ~ 0.63 are the critical exponents of order parameter 
and correlation length in the (three-dimensional) Ising universality class, and he 
related the amplitude prefactor S to an amplitude characterizing the asymmetry 
of the liquid-vapor phase diagram near criticality. However, with the methods 
described in the present paper a test of Eg. [64| seems prohibitively difficult. 

5 Wall-attached droplets and methods for esti- 
mating line tensions 

We now study the order parameter distribution of systems such as considered 
in the previous section (see e.g. Fig. [9]) but we remove the periodic boundary 
condition in z- direction such that we simulate a, L x L x D geometry with 
periodic boundaries in x- and y-directions, while in the z-direction the system is 
confined by walls under incomplete wetting conditions. If we restrict attention 
to strictly symmetric systems (Ising model, or symmetric binary LJ model), it 
is straightforward to choose a strictly "antisymmetric" boundary condition (as 
used already in Sec. 3 to obtain information on contact angles), which has the 
advantage that phase coexistence in this thin film geometry is not shifted relative 
to the bulk. Then the analysis outlined in Fig. 12 can be straightforwardly 
extended to study the free energy excess due to wall-attached droplets or slab 



configurations confined by planar interfaces (Fig. 16). For the Ising system 
shown here, the free energy excess is defined in analogy to Eq. 49 as (note that 
no bulk field H is applied here) 

)], (65) 

and the field HL{m,Hi,T) simply is the derivative of this function, in full 



analogy to Eq. 55 



HlM^,H,,T) = [dfL{m,H,,T)/dm]T,H,- (66) 

Due to the choice of the antisymmetric boundary fields, we still have a strict 
symmetry of the effective free energy and an antisymmetry of its derivative, 

fLMrn.Hi,T) = fLM-rn,H,,T), (67) 

HlM-^,H,,T) = -HlM^^H,,T). (68) 

Again, the simplest case is the slab configuration in Fig. |16[ in particular for 
the case Hi = 0, where the contact angle is 9 = tt/2. Due to the symmetry of 
our model (see Sec. 3) we have two domain walls of area LD (in d — 3), each 
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Figure 18: Surface free energy F{R,9) /ksT of wall attached droplets in the 
Ising model at ksT/J = 3.0 plotted vs. the droplet radius R, for Hi/ J = 
(0 — 90°), case (a), and Hi/ J — 0.73 {6 w 23°), case (b). Upper curve in each 
case shows 'iTTR^j{R)f{9) with 7(i?) taken from the corresponding calculation 
in the bulk. Lower curve show data obtained as L^Af {Eq. 72} , using R from 
Eq. W^ Insert shows that the difference between both curves varies linearly with 



R, hence allowing to extract a meaningful estimate of t{9), using Eq. 72 
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2 4 6 8 10 12 14 

Figure 19: Surface free energy Fs{R,9) /k^T of wall-attached droplets in the 
symmetric binary LJ mixture at T = 1, p* = 1 (crosses) plotted vs. the droplet 
radius R, for three choices of the contact angle (obtained as described in Sec. 3). 
The full curves are the corresponding predictions Fs{R,0) — A:-KR?^AB{R)f{S). 
From pM] . 
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oriented perpendicular to the walls. Thus we expect that in this case the height 
of the free energy hump in Fig. 16 can be written as 



/i,B(TO«0,i/i =0,r) = -(7 + 2r(0- -)/!?), L^<x,D 



(69) 



Taking data at fixed (large) L and varying D, this relation was tested for both 
the Ising model and the binary Lennard-Jones mixture. It was found that in 
both cases Eq. |69] is compatible with the simulation, implying that the line 
tension t{6 = 7r/2) is negative for both models. As a caveat, we mention that 



Eq. 69 disregards the problem that for the equivalent L D geometry with peri- 



odic boundary conditions instead of walls we also expect a finite size correction, 
similarly to what is seen in the geometry (Fig. |ll[ ). However, it is reasonable 
to assume that for a system with L D this "intrinsic" size effect is even 
smaller than for the geometry, and since the slope in Fig. 11 is an order of 
magnitude smaller than the slope in Fig. |17[ 3), this problem is ignored hence- 
forth. Another issue (see also the discussion of Schimmele et al. [IS]) is the fact 
that in the separation of what is attributed to be a "surface effect" and what is 
a "line effect" , certain conventions about the meaning of lengths characterizing 
the systems are inevitable. We have taken the Ising model surface area to be 
identical to the number of spins it contains, LD. However, if one would measure 
the distance from the row n = 1 to the row at n = _D as the relevant length, 
one could say the surface area is only L{D — 1). 

If one assumes that the droplets shown schematically in Fig. [T6] have a sphere- 
cap shape, the analysis based on Eqs. 58p0 {cf. also Figs. T2j T3|- can easily be 



generalized: The volume of the droplet is no longer as used in Eq. 

but reduced by a factor f{9) 
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sphcrc-cap 



= {ATrR^/3)f{6), f{9) = (2 + cos6l)(l - cos6l)V4, (70) 



and hence instead of Eq. 59 we now have (for the Ising model) 
Am = to' - m = (to" + m' )Ati f [9) / {iL^ D) 



(71) 



The same reduction factor f{9) applies to the relation replacing Eq. 60 



L^DAf = Fs{R, 9) = AnR^-fiR) + 2tiRt{9) sin 6* 



(72) 



interpreting the (absolute) excess free energy L^DAf as a total surface excess 
free energy of the sphere cap. We also allow for a correction proportional to the 
line tension t{9). Note that the factor 27r_Rsin0 is just the length of the 3-phase 
contact line of the sphere-cap shaped droplet. 



Eq. 72 makes the explicit assumption that the same function 7(-R) or ^ab{R) 
for spherical droplets in the bulk (see previous section) also applies for sphere- 
cap shaped droplets attached to a flat wall. The simulation results for both 
the Ising model and the symmetrical binary Lennard-Jones mixture in fact are 



compatible with this assumption (Figs. 18 19). The resulting estimates t{9) for 
both models are shown in Fig. [20] Note that in the Ising case (with Jg — J) one 
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t{9) plotted vs. 9, for the 
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and X in (b)t^kre derived 



Ising model with Jg/ J = 1 
mixture at T = = 1. 
from a slab geometry with 



has a second order wetting transition, and hence one expects t{9 — >• 0) — )■ [H]. 
In contrast, the LJ mixture exhibits a first order wetting transition (see Sec. 3), 
and in this case one expects t{6 — > 0) — > +oo. However, since r is negative for 
large 0, it is impUed that r must change sign before the wetting transition is 
reached [HI US]- Indeed the numerical data (Fig. 20b) are compatible with this 
expectation. 

6 Outlook 

In the present article, we have reviewed a selection of methods which have been 
used to simulate phase coexistence between fluid phases, for simple models of 
statistical mechanics, and we have discussed the information on interfacial free 
energies and related properties (contact angles at walls, line tension) that one 
can extract from these studies. Neither the detailed interfacial structure nor 
the study of wetting transitions has been discussed. However, even with this 
restricted scope we have by no means attempted to present an exhaustive cov- 
erage of the subject, but rather we have described the techniques using a few 
examples, taken from the research groups of the authors for the sake of sim- 
plicity. It needs to be stated, that very valuable research on related topics can 
be found in the literature from other groups, sometimes also for other models 
than the models (lattice gas, Lennard- Jones fluid, symmetric binary LJ mix- 
ture) that were discussed here. E.g., a particularly well-studied model system 
is the square-well fluid. Van Swol and Henderson [1541 1155j have investigated 
wetting and drying at fluid-wall interfaces for this model, and compared their 
simulations to density functional calculations. Due to the lack of any particu- 
larly symmetries for this model (or the Lennard- Jones fluid considered in the 
pioneering study by Sikkenk et al. |156| ). it was difficult to obtain very accurate 
results on interfacial free energies, contact angles etc. at this time, consider- 
ing also the fact that now orders of magnitude more computer power can be 
invested than available for this early work |1541 1155[ I156j . Significant progress 
has also been made via a better understanding of the statistical mechanics of 
the fluid- wall interactions |157j , and this work is also basic for a very interesting 
recent simulation study |158j where a slab geometry of a liquid bridge between 
two walls with contact angles 9i , 62 such that 9i + 62 ^ tt has been analyzed. 
In this work, it became possible to derive estimates for the line tension t{9) for 
several values of the contact angle 9. 

We also add that we have only addressed the problem of phase equilibria of 
two fluid phases in the presence of a solid flat wall. We have neither considered 
other geometries (e.g. wedge geometries, where "filling transitions" can occur 
|24j ) nor have we considered interfacial phenomena where the solid wall is re- 
placed by another fluid phase; in this case three fluid-fluid interfaces also can 
meet at a contact line, but the corresponding line tension |Tn| was out of con- 
sideration here. Also the problem of contact angles and wall-attached droplets 
for non- volatile fluids was out of our scope. E.g., a typical example is a droplet 
formed by a polymer melt at a substrate |160) : the vapor surrounding the poly- 



44 



mer droplet typically docs not contain any polymers at all, so strictly speaking 
the system is not in full equilibrium (which would imply that the polymers com- 
pletely evaporate), but in practice this does not occur at all at the time scales 
of interest. Some of the methods described here could be carried over to a study 
of such morc^ complex systems as well. 
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